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Preface 


This  paper  is  the  final  report  for  and  describes  the  work  accomplished  under  Grant 
No.  N00014-90-J-1294  from  the  Fluid  Dynamics  Program,  Mechanics  Division, 
Office  of  Naval  Research.  The  purpose  of  the  grant  was  to  support  research  on  a 
study  of  the  surface  layers  at  an  air-water  interface.  The  period  of  the  research  was 
15  November  1989  through  15  February  1991. 

There  are  no  other  publications  completed  tmder  this  grant;  however,  this  paper 
is  being  submitted  for  refereed  pubhcation. 
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1  Introduction 


This  paper  describes  the  application  of  a  two-equation  turbulence  model  to  a  strat¬ 
ified  two-phase  channel  flow  as  part  of  a  research  project  into  wave  generation  by 
wind.  The  project  was  prompted  by  a  review  of  wave-generation  hterature  whose 
purpose  was  to  determine  the  extent  of  the  agreement  between  extant  theory  and 
observations  of  wave  growth. 

The  hterature  review  revealed  that  wind-generated  waves  may  broadly  be  di¬ 
vided  on  the  btisis  of  u./c,  i.e.,  the  ratio  of  the  friction  velocity,  u,,  to  the  wave 
celerity,  c.  Waves  with  values  of  u,/c  >  0.15  may  be  termed  “short  waves”,  suid 
waves  with  u,/c  <  0.15  may  be  termed  “long  waves”.  The  short  wave  classifica¬ 
tion  includes  gravity-capillary  waves  and  most  wind- waves  generated  in  laboratory 
facihties.  The  long  wave  classification  includes  most  wind-generated  waves  in  the 
ocean,  as  well  as  artificially  generated  progressive  waves  over  which  wind  is  blown 
in  laboratory  facihties. 

It  was  found  that  the  growth  of  short  waves  is  quite  weU  described  by  existing 
wave  generation  theories  and  numerical  models.  For  example,  Valenzuela  ( 1976)  and 
Kawai  ( 1979)  found  very  good  agreement  between  their  experimental  measurements 
of  the  growth  of  gravity-capiUary  waves  and  the  predictions  of  viscous  hnear  stabihty 
theory. 

However,  the  growth  of  long  waves  is  only  poorly  described  by  extant  theory 
and  numericzd  models.  For  example,  wave  growth  rates  predicted  by  Miles’  (1957, 
1959)  inviscid  hnear  instabihty  theory  are  generally  2  to  5  times  smaller  than  the 
most  reliable  experimental  measurements.  Various  modifications  to  Miles’  original 
theory  have  improved  the  agreement  between  theory  and  experiment  somewhat.  For 
example,  Al-Zanaidi  eind  Hui  (1984)  found  that  the  inclusion  of  the  wave-induced 
Reynolds  stresses  in  their  numerical  model  increased  the  predicted  wave  growth 
rates  by  up  to  50%  in  some  cases.  However,  the  agreement  between  theoretical 
predictions  and  experimental  measurements  is  stiU  poor. 

Interestingly,  there  is  very  good  agreement  between  Miles’  inviscid  hnear  stabil¬ 
ity  theory  and  the  results  of  experimental  and  numerical  studies  of  the  flow  over 
sohd,  but  progressive,  monochromatic  waves  (see,  for  example,  Zagustin  et  al.,  1968 
and  Norris  and  Reynolds,  1975).  The  most  important  physical  diiference  between 
the  flow  over  an  artificial  waveform  and  the  flow  over  wind-generated  water  waves 
is  that  there  is  motion  in  the  water  under  an  actual  wind  wave,  whereas  an  ar¬ 
tificial  waveform  has  a  sohd  boundary.  This  observation  suggests  that  the  shear 
flow  in  the  water  is  important  for  the  growth  of  long  waves.  Other  evidence  also 
points  to  the  importance  of  the  shear  flow  in  the  water  for  long  wave  growth.  For 
example,  Valenzuela  (1976)  demonstrated  the  importance  of  the  shear  flow  in  the 
water  for  gravity-capillwy  (short)  waves.  He  solved  the  coupled  air- water  prob¬ 
lem  (neglecting  turbulent  effects),  and  found  that  the  inclusion  of  the  shear  flow 
in  the  water  produced  a  significant  increase  in  computed  growth  rates,  especiaUy 
for  the  waves  with  the  lowest  values  of  u,/c  (i.e.,  waves  closest  to  the  long  wave 
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classification).  Gent  and  Taylor  (1977)  found  that  the  inclusion  of  the  surface  drift 
in  their  numericaJ  model  of  wave  growth  significantly  altered  the  surface  pressure 
distribution,  but  only  produced  a  slight  increase  in  the  wave  growth  rate.  However, 
they  only  considered  two  different  cases,  and  so  the  question  remains  as  to  whether 
the  wave  growth  rate  could  change  significantly  if  a  larger  range  of  wave  param¬ 
eters  was  examined.  Thus,  it  seems  likely  that  the  inclusion  of  the  shear  flow  in 
the  water  in  a  numerical  model  of  the  growth  of  long  waves  will  increase  the  wave 
growth  rate,  which  would  improve  the  agreement  between  theoretical  predictions 
and  experimental  measurements  of  wave  growth. 

A  research  program  is  currently  imderway  to  develop  a  comprehensive  numericeJ 
model  of  wave  growth.  The  aim  of  the  model  is  to  adequately  predict  the  growth 
rates  of  both  short  amd  long  waves.  The  model  will  solve  the  coupled  tiir- water  wave 
growth  problem,  and  will  consider  the  shear  flow  in  the  water  as  well  as  including 
turbulence  effects.  As  such,  the  model  will  be  am  extension  of  previous  numerical 
models  of  wave  growth. 

The  numerical  model  is  of  the  line2ir  stability  type,  similar  to  most  models  of 
wave  growth  (e.g..  Miles,  Al-Zzmaidi  and  Hui,  and  Valenzuela).  Such  models  solve 
for  perturbations  to  a  given  meain  flow  profile  caused  by  a  single  wave  component. 
The  resulting  perturbation  pressure  and  velocity  components  may  be  used  to  com¬ 
pute  the  wave  growth  rate  and  other  quantities  of  interest. 

The  numerical  model  is  developed  in  two  parts.  The  first  involves  developing  a 
solution  for  the  equations  governing  the  mean  flow,  and  the  second  part  involves 
solving  the  perturbation  problem.  While  these  two  parts  are  essentially  separate 
problems,  it  is  important  to  consider  the  solution  of  the  mean  flow  part  in  the 
context  of  the  overall  model.  This  paper  describes  only  the  successful  solution  of 
the  mean  flow  problem. 

Many  previous  models  of  wave  growth  simply  specify  the  metin  velocity  profile 
antdytic^ily.  For  example.  Miles  (1957,  1959)  considered  a  logarithmic  velocity 
profile  in  the  air  and  ignored  the  flow  in  the  water,  whereas  Valenzuela  ( 1976)  used 
a  more  realistic  Unear-logarithmic  profile  in  the  air  and  the  water.  However,  the 
inclusion  of  turbulence  effects  in  our  model  requires  that  mean  profiles  of  turbulent 
quantities,  as  weU  as  the  velocity,  must  be  provided  as  input  to  the  perturbation 
problem.  Thus,  the  analytical  specification  of  the  mean  profiles  is  out  of  the  question 
for  our  problem,  and  it  is  essential  to  solve  the  equations  governing  the  turbulence 
and  the  mean  flow  numerically. 

It  is  weU  known  that  the  mean  flow  profile  C2in  have  a  large  effect  on  the  com¬ 
puted  growth  rate.  For  example,  van  Gastel  et  al.  (1985)  found  that  the  shape 
of  the  mean  velocity  profile  can  change  the  growth  rate  by  a  factor  of  more  than 
three.  In  addition,  the  perturbation  equations  are  sensitive  to  the  details  of  the 
turbulence  profiles  in  the  vicinity  of  the  air-water  interface.  Given  the  sensitivity 
of  the  perturbation  problem  to  the  mean  profiles,  it  is  essential  that  the  mean  flow 
equations  are  solved  through  both  boimdary  layers  at  the  air-water  interface. 

The  solution  of  the  mean  flow  problem  is  developed  for  the  case  of  a  stratified 
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two-phase  (air-water)  flow  in  the  the  Stanford  Wind  Water  W'ave  Research  Facility 
(SWWW’RF)  for  several  reasons.  First,  this  will  allow  comparisons  of  the  mean 
flow  profiles  with  previous  experimental  measurements  in  the  SWWWRF.  Second, 
this  will  also  allow  compaurisons  of  both  mean  and  wave-induced  quantities  with  the 
results  of  our  ongoing  particle  tracking  experiments  in  the  SWWWRF.  Third,  other 
test  cases  of  interest  (e.g.,  fully  developed  boimdary  layer  with  no  pressure  gradient ) 
are  easily  computed  as  simplifications  of  the  two-phase  channel  flow  problem. 

As  it  turns  out,  the  mean  flow  equations  are  quite  challenging  to  solve  numeri¬ 
cally,  since  they  are  non-hnear  aind  involve  very  steep  gradients  near  the  boundaries 
and  the  interface.  Indeed,  we  speculate  that  the  difficulty  of  solving  the  coupled 
mean  flow  equations  for  the  “long  wave”  domain  has  hindered  previous  efforts  to 
extend  the  success  of  the  “short  wave”  models  of  V'^alenzuela,  Kawai,  and  van  Gastel 
et  aJ.  However,  we  are  hopeful  that  the  considerable  effort  involved  in  solving  these 
equations  will  pay  dividends  in  the  long  rvm,  and  that  the  solution  of  the  perturba¬ 
tion  equations  using  the  computed  mean  profiles  will  yield  growth  rates  in  better 
agreement  with  experiments. 

The  following  sections  discuss  details  of  the  formffiation  of  the  model,  the  nu¬ 
merical  solution  techniques,  and  results  obtained  from  comparisons  of  the  mean 
flow  model  to  several  test  cases. 


2  Mathematical  Formulation 

2.1  Governing  equations  and  turbulence  model 

The  numerical  model  solves  for  the  turbulent,  stratified  air-water  flow  in  a  channel. 
The  flow  field  is  perturbed  by  a  two-dimensional,  monochromatic,  small  zimplitude 
gravity  wave  at  the  interface.  In  addition,  depending  on  the  windspeed,  the  inter¬ 
face  may  be  ruffied  by  short  wavelength  wind-generated  ripples.  These  ripples  are 
treated  as  roughness  at  the  interface.  The  model  is  based  on  the  Continuity  equa¬ 
tion,  the  Reynolds  averaged  Navier-Stokes  equations,  and  two  transport  equations 
for  the  turbulent  kinetic  energy  and  the  dissipation  which  are  necessary  to  provide 
closure  for  the  unknown  Reynolds  stress  terms  which  appear  in  the  Navier-Stokes 
equations.  These  equations  are  transformed  to  the  wavy  coordinate  system  of  Hsu  et 
al.  (1981)  so  that  the  water  surface  becomes  a  straight  line.  The  tr2insform  greatly 
facilitates  the  application  of  the  coupling  conditions  at  the  air-water  interface.  In 
fact,  Norris  and  Reynolds  (1975)  foimd  that  such  a  transform  was  essential  to  ob¬ 
tain  good  agreement  between  their  experiments  and  their  numerical  model.  The 
transformed  equations  are  then  expanded  in  terms  of  the  wave  slope,  ak,  where  a 
is  the  wave  amplitude,  k  =  27r/£  is  the  wave  niunber,  and  L  is  the  wave  length. 
Terms  of  order  (ait)°  in  the  resulting  equations  describe  the  mean  flow  over  the  wavy 
interface,  while  terms  of  order  (ait)*  describe  the  first  order  perturbations  to  the 
mean  flow.  As  it  turns  out,  the  zeroth  order  mean  flow  equations  in  the  transformed 
coordinate  system  are  identiced  to  the  governing  equations  for  stratified  two- phase 
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flow  in  a  channel  with  no  gravity  wave  present  on  the  interface.  Thus,  the  details  of 
the  coordinate  transformation  axe  not  presented  here.  It  suflBces  to  point  out  that 
the  solution  is  valid  either  in  an  untransformed  (cartesian)  coordinate  system  with 
no  wave  at  the  interface,  or  in  the  transformed  coordinate  system. 

The  model  assumes  that  both  fluids  are  incompressible  with  constamt  physical 
properties  and  that  the  flow  is  steady,  turbulent,  and  fully  developed.  In  our  ex¬ 
perimental  facility  the  mean  flow  does  change  slightly  with  fetch.  However,  the 
assumption  of  fully  developed  flow  is  justified  in  our  case  since  in  reality  the  mean 
flow  profile  changes  only  very  slightly  over  one  wavelength  of  the  gravity  wave  at 
the  interface.  The  well  known  k-epsilon  turbulence  model  is  used  to  model  the 
turbulence  in  the  flow.  Since  the  perturbation  equations  have  been  found  to  be 
very  sensitive  to  the  mean  velocity  profile,  it  is  essential  that  the  mean  profile  be 
computed  through  the  boimdary  layers  at  the  air- water  interface,  rather  than  using 
empirical  wzdl  functions  applied  at  some  distance  above  the  interface.  The  standard 
k-epsilon  is  only  valid  for  regions  where  molecular  viscosity  is  unimportemt  (i.e.,  fax 
from  the  boundaries  and  the  interface).  Thus,  we  have  used  an  extended  version 
of  the  standard  k-epsilon  model  that  is  valid  for  low  Reynolds  number  regions.  VVe 
decided  to  use  the  model  proposed  by  Launder  and  Sharma  (1974)  since  Patel  et 
al.  (1984)  found  that  it  performed  among  the  best  of  the  models  they  e\’aluated 
in  a  variety  of  test  cases.  This  model  has  also  been  successfully  applied  to  similar 
stratified  two-phase  flow  problems  by  Akai  et  al  (1981)  and  Issa  (1988). 

The  equations  governing  the  mean  flow  in  both  fluids  for  steady,  fully  developed 
flow  are  (Patel  et  al.,  1985) 


e  =  €  +  d 


(5) 


(6) 


k  = 


(7) 
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and 


bu'.  du\ 
^  =  1/  ~  ^ 

dxi  dxi 


(S) 


In  the  above  equations,  y  is  the  vertical  coordinate  measured  upwards  from  the 
channel  bottom,  x  is  the  streamwise  coordinate,  u  is  the  streamwise  velocity,  k  is 
the  turbulent  kinetic  energy,  and  e  is  the  (homogeneous)  dissipation  rate.  Also, 
1/  represents  the  kinematic  viscosity,  and  p  represents  the  density.  The  remaining 
constants  and  fimctions  depend  on  the  type  of  low  tmbulence  Reynolds  number 
model  used.  These  constants  and  functions  are  defined  below  for  the  Launder- 
Sharma  model; 


=  exp 


-3.4 

.(1 -I- i^^/50)^ 


(9) 

(10) 

(11) 


/.  =  1 


(12) 


/j  =  1  -  0.3 exp  (13) 

A  big  advantage  of  the  Launder-Sharma  low-Reynolds-number  model  presented 
above  is  that  the  turbulence  Reynolds  number,  Rt,  is  the  only  parameter  that 
determines  the  behavior  of  the  model  near  a  boundary.  As  a  consequence,  it  is 
possible  to  have  a  non-zero  value  of  turbulent  kinetic  energy  at  a  boimdary  with 
the  Laimder-Sharma  model.  In  contrast,  many  of  the  other  available  low  Rj  models 
were  originally  developed  for  flow  near  solid  walls,  and  require  that  the  turbulent 
kinetic  energy  is  zero  at  the  boundaries.  While  this  restriction  is  true  for  a  solid, 
smooth  wall,  it  is  not  necessarily  valid  for  rough  boundaries,  or  the  air-water  inter¬ 
face. 

For  high  turbulence  Reynolds  numbers,  Rj,  the  fimctions  and  /j  tend  to 
unity,  and  so  the  behavior  of  the  model  depends  only  on  the  values  of  the  five 
remaining  constants.  The  values  of  these  constants  are  taken  to  be  =  0.09, 
cij  =  1.44,  cjt  =  1.92,  Oft  =  1.0,  and  =  1.3,  which  are  the  same  as  for  the 
standard  high-Reynolds-number  k-epsilon  turbulence  model.  Thus,  the  Launder- 
Sharma  model  reduces  to  the  standard  k-epsilon  model  away  from  the  boundaries 
and  the  interface. 
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2.2  Application  of  the  governing  equations  to  two-phase 
stratified  channel  flow 

In  this  section  the  governing  equations  presented  above  are  applied  to  the  coupled, 
two-phase,  stratified  channel  flow  problem.  Let  H  be  the  total  height  of  the  channel, 
cind  define  the  lower  fluid  and  the  upper  fluid  to  be  layers  1  aind  2  respectively.  Let 
Q  be  the  fraction  of  the  channel  height  occupied  by  fluid  1.  Define  the  following 
non-dimensional  variables 


U  =  u/u.,  K  =  k/ul  E  =  €H/uI 
X  =  x/H,  Y  =  y/H,  P=p/(>2ul 


where  u,  is  the  friction  velocity  in  layer  2  at  the  interface.  Substituting  the  above 
non-dimensional  variables  into  the  governing  equations  results  in  the  following  equa¬ 
tions 


where 


if  0  <  r  <  o 

if  a  <  K  <  1 

if  0  <  r  <  a 
if  a  <  y  <  1 


ifO<r<a 

ifa<r<l 


=  E  =  E  +  D 


(IS) 

(19) 


(20) 

(21) 


The  quantity  E  =  iH/ul  is  the  “dissipation  variable”.  We  solve  a  tramsport 
equation  for  E  (rather  than  E)  for  numerical  convenience,  since  E  is  a,  known 
quantity  at  a  boimdary  (see  Patel  et  al.,  1984). 

The  sixth  order  set  of  ordinary  differential  equations  presented  above  requires 
six  boimdary  conditions.  The  equations  also  require  two  auxiliary  conditions  in 
order  to  select  the  pressure  gradients  in  each  layer,  as  well  as  appropriate  coupUng 
conditions  at  the  interface  of  layer  1  and  2. 

The  upper  and  lower  walls  may  be  considered  as  hydraulically  smooth  in  our 
facihty,  so  that  the  boimdary  conditions  at  the  lower  and  upper  walls  are  given  by 
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K  =  0,  £  =  0,  U  =  0  at  K  =  0  and  Y  =  I  ( 23 ) 

Coupling  conditions  are  required  to  relate  the  variables  in  layer  1  to  those  in 
layer  2  at  the  interface.  For  our  problem,  it  is  known  that  both  the  velocity  and 
stress  cire  continuous  across  the  interface.  However,  similar  relationships  are  not 
readily  apparent  for  the  turbulent  kinetic  energy  and  the  dissipation  rate  (Issa. 
1988).  Thus,  it  was  decided  to  solve  a  quasi-coupled  problem  in  which  values  of  K 
and  E  are  prescribed  at  the  interface  in  a  similar  way  to  a  sohd  boundaiy,  whereas 
the  velocity  and  the  stress  are  assumed  to  be  continuous  there.  The  stress  and 
velocity  couphng  conditions  at  F  =  a  are  given  by 
velocity  continuity; 

l\=U2  (24) 

stress  continuity: 


The  prescribed  values  of  K  and  E  may  be  expressed  as 


For  a  smooth  interface,  A'o^,  Eq^,  and  Eqw  zero,  whereas  for  a  rough 

interface  these  values  aure  non-zero.  We  found  that  the  specification  of  a  non-zero 
value  of  turbulent  kinetic  energy  at  the  interfaice  as  used  by  Akai  et  ad.  (1981) 
cuid  Issa  ( 1988)  could  not  predict  the  measured  velocity  profiles  in  our  channel  for 
free  streaun  wind  speeds  greater  than  about  3  m/s.  Thus,  we  developed  a  more 
rational  method  for  predicting  the  boundary  vadues  of  turbulent  kinetic  energy*  amd 
dissipation  rate  at  a  rough,  wavy  interface  that  is  appUcable  to  smooth,  trainsitionad. 
and  fully  rough  flow  regimes.  It  is  assumed  that  the  interfaciad  roughness  may  be 
chairaw:terized  by  a  roughness  length,  yo,  which  is  defined  as  the  vertical  intercept  of 
a  semilogairithmic  plot  of  y  versus  u(y)  —  u,,  where  u,  is  the  surface  drift  velocity 
velocity  (Amorocho  amd  DeVries,  1980).  Both  the  water  and  the  air  are  assumed 
to  be  subject  to  the  same  roughness  elements  (Kondo,  1976),  so  that  the  roughness 
Reynolds  number  in  the  adr  is  about  twice  that  in  the  water.  Effectively,  this  means 
that  the  interfaiciad  region  is  treated  in  a  similar  manner  to  a  sohd  boundary.  That 
is,  all  effects  of  the  orbital  motions  of  the  water  waves  have  been  ignored  for  the 
present  time. 

The  vadues  of  K  and  E  at  the  interface  may  be  determined  by  ensuring  that  the 
eddy  viscosity  at  the  wall  (given  by  Equation  [21])  is  consistent  with  that  predicted 
by  the  mixing  length  model,  ais  well  as  by  a  one-equation  turbulence  model.  The 
eddy  viscosity  for  the  mixing  length  model  is  given  (in  non-dimensionad  form)  by 
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(Kays  and  Crawford,  1980),  whereas  the  eddy  viscosity  for  a  one-equation  turbu¬ 
lence  model  is  given  by 

(28) 

(ASCE  Task  Committee  on  Turbulence  Models  in  Hydraulic  Computations,  1988) 
where  L  is  the  turbulent  length  scale  (or  mixing  length).  By  noting  that  Lq  =  nyo/H 
at  the  interface  (K  =  a)  and  eliminating  ut  from  Equations  [21],  [27],  and  [28].  the 
following  expressions  for  Kq  and  are  obtained 

ifo  =  (-v/fl  +  ^(VIRf-ViLlIv)  (29) 

=  *)/!«.  (30) 

Thus,  if  the  non-dimensional  roughness  height,  Yq  =  yo/H,  is  known,  the  values  of 
Kq  and  Eq  on  either  side  of  the  interface  (i.e.,  ^a) 

determined.  In  practice,  the  roughness  height,  yo,  may  be  related  to  the  free  stream 
wind  speed  and  the  friction  velocity,  u.  (Amorocho  aind  DeVries,  1980). 

The  values  of  Kq  and  Eq  predicted  from  Equations  [29]  and  [30]  have  the  correct 
assymptotic  behavior.  That  is,  for  Iq  — 0  (i.e.,  smooth  boimdary)  they  yield 
A'o  =  0  and  Eq  =  0  at  the  interface,  whereas  for  large  Yq  (i.e.,  fully  rough  flow) 
they  predict  Kq  — ♦  l/(Py/^)  and  Eo  — ►  l/(«lo^^*),  which  are  exactly  the  wall 
functions  applied  at  K  =  Iq  (ASCE  Task  Committee  on  Turbulence  Models  in 
Hydraulic  Computations,  1988).  Thus,  they  provide  a  general  method  for  specifying 
boundary  conditions  on  the  turbulent  kinetic  energy  and  the  dissipation  rate  at  a 
smooth,  transitional,  or  rough  boundary. 

Finally,  conditions  are  required  to  set  the  pressure  graidients  in  each  layer.  For 
fully  developed  channel  flow  the  pressure  gradient  in  each  layer  must  be  a  constant 
value.  However,  these  pressure  gradients  are  initially  unknown.  Also,  since  the  flow 
in  each  layer  is  coupled,  the  pressure  gradient  in  layer  1  affects  the  flow  in  layer  2 
and  vice  versa.  Since  there  is  no  net  flow  in  the  lower  layer  for  a  laboratory  facility 
such  as  the  SWWWRF,  the  pressure  gradients  must  be  chosen  such  that  the  average 
velocity  (i.e.,  flow  rate)  in  layer  1  is  zero.  Another  condition  that  must  be  satisfied 
by  the  solution  is  that  the  non-dimensional  stress  at  the  interface  in  layer  2  must 
be  unity  (this  is  a  consequence  of  normalizing  the  solution  on  the  friction  velocity 
at  the  interface  in  layer  2).  These  conditions  may  be  expressed  as 

rUiY)dY  =  0  (31) 

Jo 

Given  the  above  condition  on  the  stress  at  the  interface.  Equation  [15]  may  be 
integrated  to  yield 

f{Y)  =  -  a)  +  1  (33) 
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where 


(34) 


is  the  total  shear  stress  normalized  by  This  integration  reduces  the  order 

of  the  equation  set  by  one.  Since  the  auxiliary  condition  relating  to  the  stress  at 
the  interface  has  been  used,  the  new  auxiliary  condition  on  the  choice  of  pressure 
gradient  is  chosen  to  be  the  boundary  condition  17(1)  =  0. 

3  Numerical  solution 

The  governing  equations  to  be  solved  are  a  5th  order  set  of  non-linear  ordinary  differ¬ 
ential  equations.  Along  with  the  boundary  and  coupling  conditions,  they  comprise 
a  3  point  boundary  value  problem  (BVP).  In  our  case  the  problem  is  a  difficult  one. 
since  it  involves  the  computation  of  profiles  through  four  boundary  layers  -  one  at 
each  wzdl,  and  two  at  the  air-water  interface.  In  addition,  the  low-Reynolds-number 
k-epsilon  model  introduces  increased  computational  difficulty. 

Initially,  we  attempted  to  solve  the  problem  using  the  shooting  method.  This 
method  had  been  successfully  applied  to  a  simpler  channel  flow  problem  by  Nor¬ 
ris  and  Reynolds  (1975),  who  used  a  one-equation  turbulence  model.  Briefly,  this 
method  involves  guessing  unknown  boundary  conditions  at  the  channel  bottom, 
then  integrating  (or  shooting)  toward  the  top  wall  where  the  known  boimdary  con¬ 
ditions  are  checked.  If  the  top  wall  boundary  conditions  are  not  satisfied,  then 
a  new  guess  is  made  for  the  bottom  boimdary  values,  and  the  process  continues 
until  the  top  wall  boundary  conditions  are  satisfied.  Unfortimately,  the  governing 
equations  proved  extremely  sensitive  to  small  changes  in  the  two  unknown  bottom 
boundary  conditions  (i.e.,  the  equations  are  very  stiff  in  the  mathematical  sense), 
md  so  the  shooting  method  failed  for  our  problem.  On  the  other  hand,  Norris  and 
Reynolds  only  had  one  unknown  boundary  condition  at  the  channel  bottom;  thus, 
it  was  possible  to  use  the  shooting  method  effectively  in  their  case,  even  though 
their  problem  was  still  very  sensitive  to  the  unknown  boundary  condition. 

Next,  the  equations  were  solved  by  the  method  of  Patankar  and  Spalding  (1967) 
for  parabolic  partial  differential  equations.  This  method  involved  reintroducing 
the  time  derivatives  into  the  governing  equations,  finite  differencing  the  resulting 
partial  differential  equations,  then  marching  forward  in  time  from  an  initial  guess 
to  steady  state.  The  time  meu'ching  technique  chosen  is  second-order  accurate  and 
implicit,  so  that  the  method  is  imconditionally  stable.  In  our  case  the  non-linear 
equations  must  be  linearized  at  each  step.  In  practice,  this  linearization  limits  the 
time  step  that  can  be  used  to  achieve  convergence  to  a  small  value.  If  larger  time 
steps  are  used,  the  solution  may  become  non-physical  (e.g.,  negative  values  of  k  and 
e).  This  method  of  solution  proved  to  give  satisfactory  results.  However,  the  rate 
of  convergence  towards  steady  state  is  extremely  slow  in  some  cases,  and  a  large 
number  of  iterations  (e.g.,  10,000  to  20,000)  are  required  to  reach  steady  state.  This 
slow  rate  of  convergence  was  also  found  by  Akai  et  al.  (1981)  who  used  a  similar 
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time  maxdung  method  to  solve  the  stratified  two-phase  flow  problem,  and  found 
that  about  10000  iteration  cycles  were  required  to  reaoh  steady  state. 

Finally,  the  problem  was  solved  using  a  conceptually  simple  finite  difference 
method  for  non-hneau'  BVPs  outhned  by  Ascher  et  ad.  (1988).  With  this  technique, 
an  initial  guess  at  the  solution  is  updated  iteratively,  using  Newton  iteration,  until 
convergence  is  achieved.  This  method  has  been  found  to  give  extremely  rapid 
convergence  to  the  final  solution  (theoreticadly,  the  convergence  is  quadratic).  For 
example,  given  a  suitable  initial  guess,  the  problem  converges  to  the  solution  in 
only  10  to  20  iterations.  The  excellent  convergence  properties  of  Newton's  method 
do  not  come  without  a  price,  however.  First,  the  method  requires  the  specification 
of  the  Jacobian  matrix  which  involves  very  comphcated  algebraic  expressions  for 
our  particular  problem.  Second,  the  method  requires  a  good  initial  guess  to  the 
solution  (a  known  characteristic  of  Newton’s  method).  Fortunately,  a  good  initial 
guess  was  available  in  our  case,  since  we  had  already  solved  the  problem  using  the 
time-marching  technique.  Once  a  solution  was  obtained  for  one  test  case  (e.g.,  for 
a  given  Reynolds  number),  it  was  used  as  an  initial  guess  for  related  test  cases. 
Deteiils  of  the  application  of  the  finite  difference  method  with  Newton  iteration  to 
our  problem  are  presented  in  the  following  section. 

3.1  Reduction  of  the  problem  to  a  set  of  first  order  ODEs 

The  BVP  was  reduced  to  a  coupled  set  of  first  order  ODEs  for  solution  by  the  finite 
difference  technique.  This  reduction  has  two  main  advantages.  First,  the  gradients 
of  the  principal  variables  (i.e.,  U,  K,  and  E)  are  obtained  from  the  solution,  as  well 
as  the  principal  variables  themselves.  Since  the  perturbation  equations  require  mean 
profiles  and  their  first  and  second  derivatives,  this  means  our  numericsd  profiles  only 
have  to  be  differentiated  numerically  once,  instead  of  twice,  thus  ehminating  one 
possible  source  of  error.  Second,  the  reduction  technique  greatly  facilitates  the 
application  of  greidient  (Neumeinn)  boimdary  conditions  eind  coupling  conditions, 
which  were  used  for  some  test  cases. 


By  defining  Ac;  =  V/R  -k-  ut,  \k  =  V/R  -b  vrjdk-,  and  \e 
may  write  the  above  equations  as 

=  F/i2  -f  we 

z'  =  f(r,z)  o<r<i 

(35) 

where 

(36) 

f  =  {A/AK,5/Ae,C/Ai;,-5K,-5Ef 

(37) 

and 

SK  =  Pr(£)'-£-D 

(38) 
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E  E^ 

Sb  =  cuA^Pr  -  c„h-  +  G 

(39) 

II 

(40) 

dE 

(41) 

dU 

(42) 

The  above  equations  are  solved  using  the  distribution  o{  C{Y)  found  from  Equa¬ 
tion  [33],  along  with  boundary,  coupling,  and  auxiliary  conditions  given  by  Equa¬ 
tions  [23]  to  [32]. 

3.2  Finite  difference  method 

The  finite  difference  method  with  Newton  iteration  for  non-linear  problems  given  by 
Ascher  et  al  ( 1988)  is  presented  below,  using  their  notation.  Consider  our  non-linear 
boundatry  value  problem  (BVP)  of  order  iV  =  5  given  by 

z'  =  f(r,z)  o<r<i 

with  boundary  conditions 

g(z(0),z(l))  =0 

For  numerical  approximation  consider  a  mesh  represented  by 

TT :  0  =  E]  <  <  ...  <  Ym+i  =  q  =  1^+2  <  ...  <  Y2M+2  =  1  (•15) 

where  1  <  j  <  M  -f- 1  for  layer  1 ,  and  M  +  2  <  i  <  2M  ■¥  2  for  layer  2.  The  points 
A/ -1-1  and  M-\-2  are  the  interfacial  grid  points  in  layer  1  and  2,  respectively.  Denote 
the  vector  of  approximate  solution  values  at  the  mesh  points  by  z,. 

A  one-step  scheme  (i.e.,  a  scheme  in  which  the  difference  operator  is  baised  only 
on  values  related  to  one  subinterval  of  the  mesh  at  a  time)  is  used  to  represent 
the  derivative  term.  One-step  schemes  have  the  considerable  advantage  over  the 
more  conventional  (at  least  in  fluid  mechanics)  two-step  schemes  in  that  the  order 
of  accuracy  of  the  finite  difference  method  is  preserved  for  a  non-uniform  grid.  That 
is,  if  the  one-step  scheme  is  second  order,  the  finite  difference  equations  will  also  be 
second  order,  regardless  of  whether  the  grid  spacing  is  uniform  or  not.  Two  of  the 
most  well  known  one-step  methods  are  the  midpoint  scheme  and  the  trapezoidal 
method,  which  are  both  second  order  accurate.  Our  BVP  exhibits  a  mild  singularity 
at  the  boimdaries,  since  /if  =  0  and  E  =  0  results  in  Pr,  D  and,  hence,  f  being 
undefined.  Thus,  for  otir  problem,  the  midpoint  scheme  is  particularly  appropriate, 
since  it  does  not  actually  evaluate  the  vector  f  at  the  grid  points  -  only  in  between. 
Application  of  the  midpoint  scheme  to  the  ODE  yields 


(43) 

(44) 
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(46) 

with  boundary  conditions  given  by 

g(*l,*2A/+2)  =  0 

(47) 

where  l^+i/2  =  ^i+i  ~  Yi-  Since  there  is  a  discontinuity  in  the  fluid  properties  at 
the  interface,  the  finite  difference  equation  above  does  not  hold  for  the  grid  point 
i  =  A/  +  1.  Instead,  analytical  relationships  (i.e.,  the  coupling  conditions)  relate 
the  interfacial  variables  in  layer  1  to  those  in  layer  2. 

Equations  [46]  and  [47],  along  with  the  interfacial  coupling  conditions,  repre¬ 
sent  a  system  of  N{2M  -f  2)  non-linear  equations  for  the  N{2M  -f-  2)  unknowns, 
z,.  Newton’s  method  is  used  to  solve  the  non-linear  system.  Application  of  New¬ 
ton’s  iteration  method  to  the  above  non-linear  problem  results  in  (see  Ascher  et  al. 

(1988)) 

- 5 +  ^(K+,/2)w,]  =  -y,zT 

(48) 

where 

iV,zr  =  ~  f{Y,^v/2, +  *.>.)) 

(49) 

A{YiJrii2)  =  ^  ^^+1/2,  +*•+1)^ 

(50) 

Here  z”  is  the  known  solution  after  m  iterations  (z°  is  the  initial  guess), 
boundary  conditions  at  the  lower  ^lnd  upper  walls  may  be  written  as 

The 

S,Wi  -H  B(,W2M+2  =  -g(*r’*2M+2) 

(51) 

where 

dzT 

(52) 

and 

0  _  ^g(*ri*"w+2) 

^2M+2 

The  next  iterate  is  determined  by 

(53) 

=  y”*  -f  W. 

(54) 

The  hnear  system  represented  by  Equation  [48]  may  be  written  as 

5,w.  -H  ii,w,+,  =  q, 

(55) 

where  Si  and  Ri  are  the  N  x  N  matrices 

(56) 
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iZ, —  (57) 

and  qi  is  the  right  hand  side  vector  formed  from  the  cinrently  available  solution 
vector  Z-"  given  by 

q.-  =  (58) 

For  our  problem  the  boundary  conditions  at  the  upper  and  lower  walls  are 
separated.  Suppose  there  are  nu  boimdary  condtions  at  the  upper  boundary.  Then, 
there  must  be  nl  =  N  —  nu  boundary  conditions  at  the  lower  boundary  (for  our 
problem  nu  =  2  and  nl  =  3).  By  selecting  the  appropriate  rows  from  the  matrices 
Bo  and  Bt,  the  boundzury  conditions  may  be  expressed  as 

B/w,  = -g,(z7*)  (59) 

and 

BuW2M+2  =  ~8«*(*2M+2)  (^0) 

where  Bu  is  an  nu  x  iV  matrix,  Bj  is  an  nZ  x  iV  matrix,  is  a  vector  of  length  nu, 
and  gt  is  a  vector  of  length  nl.  Finally,  the  coupling  conditions  at  the  interface  may 
be  written  as 

•S'aWM+l  -I-  BaWw+2  =  Qa  (61) 

The  detailed  matrix  structure  of  the  linear  system  is  shown  below: 


For  compactness,  this  equation  is  expressed  as 

Aw,  =  0  (63) 

3.3  Solution  algorithm 

At  each  iteration  step,  Equation  [63]  above  must  be  solved  to  obtain  the  vector  w,. 
One  method  of  solving  the  large  staircase  matrix.  A,  is  to  consider  it  as  banded, 
with  a  lower  bandwidth  o{  N  nl  —  and  an  upper  bandwidth  of  iV  +  nu  —  1. 
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Gaussizui  elimination  requires  an  increase  in  the  upper  bandwidth  to  2N  -  1,  and  so 
storing  the  matrix  requires  an  array  of  size  N{2M  +  2)  x  {3N  +  n/  —  1).  This  means 
that  the  storage  requirements,  as  well  as  the  computational  effort  required  to  solve 
the  matrix  are  optimally  linear  in  the  number  of  grid  points,  2M  +  2.  Although 
there  are  more  efficient  methods  available  to  solved  this  staircase  matrix  (Ascher  et 
ad.,  1988,  Chapter  7),  the  banded  method  was  chosen  since  routines  to  solve  banded 
matrices  are  readily  available  in  standard  subroutine  libraries.  For  our  problem,  the 
routines  DGBFA  and  DGBSL  from  the  LINPACK  subroutine  library  were  used  to 
factor  and  solve  the  banded  matrix  at  each  step  in  the  solution  procedure. 

Given  a  Reynolds  number,  iZ,  and  two  trial  pressure  gradients,  |y,  and  Iyj.  the 
solution  to  the  boundary  value  problem  may  be  found  using  the  following  algorithm 
(Ascher  et  al.,  1988,  Chapter  5). 

1.  Generate  a  mesh,  tr,  and  an  initial  guess,  2°.  Decide  on  a  solution  tolerance,  TOL. 
REPEAT 

2.  Generate  Bu,  Bi,  g^,  and  gi. 

3.  FOR  i  =  1,...,2A/  +  1  DO 

IF  z  =  a  THEN 

Generate  Sa,  Ra,  and 

ELSE 

Generate  Si,  R,,  and  q, 

4.  Solve  ,4w,  =  S  for  w,. 

5.  FOR  i  =  1 . 2A/  +  1  DO 

UNTIL  I'w,)  <  TOL. 

For  our  problem  the  solution  tolerance  was  set  to  10“*  or  less.  This  value  is  some¬ 
what  conservative  as  the  changes  in  the  solution  were  negligible  for  values  of  |Wt| 
less  than  about  10”^.  However,  only  one  or  two  iterations  were  required  to  re¬ 
duce  the  value  of  |w,|  from  10”^  to  10~®,  since  the  convergence  of  the  method  is 
approximately  quadratic. 

The  most  difficult  part  of  the  solution  solution  procedure  for  our  problem  is 
the  determination  of  the  Jacobian'matrix,  A,  defined  by  Equation  [50],  since  some 
components  of  f  involve  the  solution  variables  z  in  a  very  complicated  manner.  The 
derivative  terms  which  comprise  the  5x5  matrix  A  were  evaluated  using  the  sym¬ 
bolic  manipulation  program  Mathematical^ ,  which  greatly  simplified  the  process. 
The  resulting  algebraic  derivative  expressions  were  output  from  Mathematical in 
FORTRAN  form  and  incorporated  directly  into  a  subroutine  so  that  translation 
errors  were  avoided. 

For  the  boundary  conditions  presented  in  Section  2.2  the  boimdary  condition 
vector,  g,  may  be  written  as 

g=  {A'(0),A'(l),E(0),B(l),U(0)f  (64) 
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ajid  the  solution  vectors  at  the  boundaries  2ure 

I,  =  {A-(0),£(0),U(0),^(0),fl(0)}’' 

Thus,  the  boundary  condition  matrices  are  given  by 

/  1  0  0  0  0  \ 

0  0  0  0  0 

5,  =  0  10  0  0 

0  0  0  0  0 

^  0  0  1  0  0  ; 

amd 

/  0  0  0  0  0  \ 

1  0  0  0  0 

Bb=  0  0  0  0  0  (68) 

0  10  0  0 

V  0  0  0  0  0  > 

from  which  the  matrices  Bi  amd  Bu,  amd  the  vectors  g;  and  aire  readily  appairent. 

The  solution  obtained  using  the  above  algorithm  will  not,  in  general,  satisfy  the 
auxihau'y  conditions  given  by  Equations  [23]  and  [31].  Thus,  the  pressure  gradi¬ 
ents  are  adjusted,  amd  the  problem  re-solved  imtil  the  two  auxiliary  conditions  are 
satisfied.  The  aidjustment  of  the  pressure  grauiients  to  meet  the  auxiliary  condi¬ 
tions  is  analogous  to  finding  the  roots  two  non-linear  equations.  This  process  is 
formalized  as  follows.  Define  the  following: 


II 

H 

mp 

(69) 

dP 

dX2-^^ 

(70) 

/i(x.,x,)=  rU{Y)dY 

Jo 

(71) 

amd 

hixx.xt)  =  U{\) 

(72) 

The  problem  is  to  determine  (ii,X2)  such  that  fx  amd  /j  are  zero. 
fx  amd  /j  may  be  expanded  in  a  Taylor  series  au 

The  functions 

/i(x  -1-  ^x)  =  /,(x)  -1-  Y, 

;=1 

(73) 

Since  we  want  /^(x  -1-  6x)  =  0,  then  we  may  write 

2 

Y  =  0i 

J=l 

(74) 

(65) 

(66) 

(67) 
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where 


"0  =  ^  andA  =  -/<  (75) 

The  solution  algorithm  for  finding  /i  and  /j  is  as  follows.  Given  an  initial  guess 
for  (xijXj),  the  BVP  is  solved  to  yield  /i  and  /j.  Then,  the  initial  guesses  are  each 
perturbed  by  a  smzdl  amount,  and  the  problem  re-solved  to  yield  the  fimctions  ( /i  -I- 
1  /a  +  corresponding  to  (xt  -f  6xi ,  Xj)  and  (/i  +  ,  /j  -I-  )  corresponding 

to  (xi,xj  -t-  Sx-i).  The  Jacobian  terms  are  estimated  as  follows: 


a/j  sf: 

0122  —  ^ ^  T -  (79) 

UX\  0X2 

The  correction  to  be  applied  to  the  pressure  gradients  may  be  determined  by  solving 
the  linear  system  given  by  Equation  (74)  for  the  vector  {6x1, 6x2).  The  process 
is  repeated  until  fi  and  /j  are  sufficiently  close  to  zero.  For  our  problem,  the 
adjustment  of  the  pressure  gradients  to  satisfy  the  auxiliary  conditions  usually  only 
required  about  1  or  2  cycles  to  reduce  the  absolute  errors  in  fi  and  /2  to  less  than 
10“^ 


4  Results 

4.1  Model  parameters  and  initial  conditions 

The  results  of  the  numerical  model  presented  in  the  previous  sections  are  compared 
to  experimental  results  for  the  case  of  air-water  flow  in  the  Stanford  Wind  Water 
Wave  Research  Facility  (SWWWRF).  The  density  and  viscosity  of  air  were  t2iken 
to  be  1.2  kg/m^  and  1.5  x  10”*  m*  s”*  respectively,  while  those  for  water  were 
taken  to  be  1000  kg/m^  and  1.0  x  10”*  m’ s”*.  These  values  result  in  the  following 
parameters  for  the  numerical  model: 

^  =  833.3  (SO) 

—  =0.0667  (81) 

i/j 

The  total  channel  depth,  H,  is  1.93  m,  and  the  water  depth  for  the  experiments  was 
0.965  m,  so  that  a  =  0.5.  During  the  experiments  the  water  stirface  tilts  slightly, 
causing  the  reverse  pressure  gradient  in  the  lower  layer,  |yj.  However,  the  change 
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in  water  surface  level  is  very  small  (of  the  order  of  a  few  mm)  and  so  a  was  set  at 
0.5  for  all  cases.  Since  p\l >  1  for  air- water,  the  stress  at  the  upper  wall  and 
the  interface  axe  approximately  equal  in  the  air  phase.  Also,  the  stress  at  the  lower 
wall  is  approximately  zero  in  the  water  phase.  This  results  in  initial  estimates  of 
the  pressure  gradients  in  the  air  and  water  as 


dXx 


(82) 


and 


dP 

dXi 


-2 


=  -4.0 


(83) 


(1-a) 

These  pressure  gradients  are  adjusted  during  program  execution  using  the  method 
described  in  Section  3.3. 


4.2  Comparison  of  model  results  with  other  studies 

As  discussed  in  the  introduction,  Valenzuela  ( 1976)  solved  the  gravity-capillary  wave 
growth  problem  for  a  coupled  shear  flow.  He  found  excellent  agreement  between 
his  numerical  model  and  experimental  investigations  of  the  growth  of  short  waves. 
One  of  the  key  factors  in  Valenzuela’s  success  was  his  use  of  a  coupled  model  which 
took  into  account  the  effect  of  the  shear  flow  in  the  water.  He  used  an  analytical 
logarithmic-linear  velocity  profile  to  represent  the  coupled,  turbulent  mean  flow 
field.  In  the  air  the  profile  is  given  by 


u(y)  = 


( 


+  Kvlvo 

u,  -f  tij  -(-  2.5u“(ao  —  tanhao/2) 


for  y<y\ 
for  y>y\ 


(84) 


emd  in  the  water  the  profile  (in  which  y  is  positive  into  the  water)  is  given  by 


u{y)  = 


{ 


Uf 


<y/yo  for  y  <  yr 

—  2.5ulJ'(a,„  —  tanhQ,„/2)  for  y  >  yj" 


(85) 


wii6r6 

sinha,  =  (y  -  y?)/1.25y5  sinha,„  =  (y  -  yr)/l-25y^ 

yg  =  ^  ’ 

In  the  above  equations,  u,  represents  the  drift  velocity  at  the  interface,  and  su¬ 
perscripts  “a”  and  “w”  refer  the  to  air  and  water,  respectively.  Valenzuela  used 
=  0.8  (i.e.,  u,/u"  =  23.1),  and  considered  two  profiles  specified  by  yi/yo  =  5 
and  yi  /yo  =  8. 

Figure  1  shows  a  comparison  of  Valenzuela’s  mean  velocity  profile  in  the  water 
with  the  numerical  model  results  for  the  case  of  a  smooth  interface  and  Uqo  =  1-5 
m/s.  The  agreement  between  the  numerical  results  and  the  analytical  profile  is 
excellent.  The  small  difference  between  the  numerical  model  and  the  analytical 
profile  near  the  interface  is  because  the  model  predicts  u,/ul  =  0.77  instead  of  0.8 
used  by  Valenzuela.  If  Valenzuela’s  model  is  corrected  for  the  drift  velocity,  it  turns 
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out  that  the  numerical  profile  falls  between  the  two  cases  chosen  by  Valenzuela, 
which  indicates  that  the  Valenzuela  profiles  were  appropriate  bounds  to  the  actual 
profile.  The  agreement  with  Valenzuela’s  profile  is  encouraging.  In  fact,  one  might 
argue  that  the  coupled  numericzd  model  is  wasted  since  the  numerical  results  are  so 
close  to  the  assumed  profile  of  Valenzuela.  However,  this  is  not  the  case.  The  reason 
is  that,  for  gravity  wave  growth,  turbulence  effects  are  important  (as  opposed  to  the 
growth  of  gravity-capilliaxy  waves  simulated  by  Valenzuela).  Thus,  the  perturbation 
equations  require  mean  profiles  of  turbulent  kinetic  energy  and  dissipation  as  input, 
in  addition  to  the  mean  velocity  profile.  Analytical  expressions  for  these  profiles 
of  mean  turbulence  quantities  are  not  available,  and  so  our  numerical  model  for 
the  mean  flow  is  essential.  In  addition,  the  analytical  velocity  profile  given  by 
Valenzuela  does  not  have  continous  second  derivatives  which  causes  a  problem  in 
the  perturbation  solution. 

The  numericzd  model  was  applied  to  the  coupled  air- water  flow  in  the  SVV'WWRF. 
cind  the  results  compared  with  previously  measured  experimental  data.  .\11  the  ex¬ 
perimental  data  were  measured  at  a  fetch  of  13  m  from  the  air  inlet.  Thus,  each 
experiment  may  be  characterized  solely  by  the  free  stream  wind  speed.  The  mea¬ 
sured  friction  velocities,  channel  height,  and  kinematic  viscosity  of  air  were  used 
to  determine  the  Reynolds  numbers  for  the  simulations.  The  roughness  height,  Vo. 
was  determined  from  experimental  wind  profiles  using  the  relationship 


ro  = 


(ST) 


H  exp  u{2)k/u, 

in  which  the  vertical  coordinate  z  is  measured  from  the  air-water  interface.  The 
above  expression  follows  directly  from  the  well  known  logarithmic  velocity  profile 
for  rough  boundaries  (Amorocho  and  DeVries,  1980). 

Figure  2  shows  a  plot  of  the  velocity  profile  in  the  air  for  a  free  stream  w'ind  speed 
of  Uoo  —  3.2  m/s,  and  Figure  3  shows  the  corresponding  velocity  profile  in  the  water. 
The  non-zero  velocity  at  the  interface  may  clearly  be  seen.  This  drift  velocity,  u,, 
was  found  to  be  approximately  3  percent  of  the  wind  speed  in  all  cases.  The  return 
flow  in  the  water  caused  by  the  reverse  pressure  gradient  is  also  evident  in  Figme 
3.  Velocity  profiles  for  other  wind  speeds  are  similar  to  the  profiles  presented. 

The  numerical  model  was  used  to  simulate  the  experimental  results  of  McIntosh 
et  al.  ( 1975)  who  measured  the  velocity  profiles  above  the  air- water  interface  in  the 
SWWWKF  using  a  pitot  static  tube.  Figure  4  shows  a  comparison  of  McIntosh’s 
data  with  the  numerical  simulation  for  wind  speeds  of  2.3  m/s,  3.3  m/s,  5.4  m/s, 
8.1  m/s,  and  10.7  m/s.  The  agreement  between  the  numerical  model  and  the 
experimental  points  is  excellent  for  the  four  lowest  wind  speeds,  but  the  model 
underpredicts  the  roughness  Reynolds  number  of  the  highest  wind  speeds  somewhat. 
However,  the  simulation  is  still  within  the  error  bound  of  the  experimental  data 
(about  ±10  percent  in  u+),  even  at  the  highest  wind  speed. 

Next,  the  model  was  used  to  simulate  the  experimental  results  of  Cheung  (1985, 
1988)  who  measured  the  velocity  profile  in  the  water  beneath  the  air-water  interface 
using  a  Laser  Doppler  Anemometer.  Figure  5  shows  a  comparison  of  Cheung’s  data 
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with  our  model  results  for  7  wind  speeds  ranging  between  1.5  m/s  and  13.1  m/s. 
The  quaditative  agreement  between  the  measured  and  predicted  profiles  is  excellent, 
showing  that  the  intercept  of  the  profiles  decreases  with  increasing  wind  velocity 
as  the  water  surface  becomes  rougher.  The  quantitative  agreement  between  the 
experimental  profiles  and  the  numerical  simulation  may  be  seen  to  be  quite  good 
for  the  three  highest  and  the  two  lowest  wind  speeds.  However,  the  model  does  not 
predict  the  3.2  m/s  or  4.7  m/s  data  as  well.  One  re^lson  for  the  apparent  quantitative 
difference  may  be  attributed  to  uncertainty  in  the  experimental  data.  There  is  a 
large  uncertainty  involved  in  measuring  the  surface  drift  velocity  which  translates 
into  a  large  imcertainty  in  the  intercept  of  the  experimental  profiles  in  the  water. 
For  example,  consider  the  3.2  m/s  case.  Cheimg  quotes  an  imcertainty  of  ±\0% 
for  u,  which  corresponds  to  8.4  mm/s.  Dividing  by  the  friction  velocity,  u,  =  4.9 
mm/s  means  that  the  intercept  of  the  3.2  m/s  profiles  may  be  shifted  up  or  down 
by  u+  =  1.7.  In  addition  to  this  uncertainty,  there  is  at  least  a  ±10%  uncertainty  in 
the  friction  velocity,  u.,  amd  a  ±2%  error  in  the  velocity  measurements  themselves. 
Thus,  the  combination  of  these  errors  results  in  a  large  uncertainty  in  u+.  and  so 
the  failure  of  the  expermental  profiles  to  quantitatively  agree  with  the  predicted 
profiles  in  all  cases  is  not  surprising. 

It  can  be  seen  that  the  slope  of  the  experimental  data  in  the  water  deviates 
from  the  expected  value  of  l//c  for  low  values  of  y+.  This  deviation  is  caused  by  the 
effects  of  wave-turbulence  interaction  (see  Howe  et  al.,  1982)  disrupting  the  usual 
balance  of  production  and  dissipation  rate  that  occurs  in  the  vicinity  of  a  solid 
boundary.  Since  our  niunerical  model  does  not  consider  this  effect,  the  deviation  of 
the  experimental  data  points  from  the  simulation  is  to  be  expected  for  low  values 
of  t/+  (i.e.,  close  to  the  interface). 

Figure  6  shows  profiles  of  turbulent  kinetic  energy  predicted  by  our  model  in 
the  air  above  the  interface  over  a  range  of  wind  speeds.  The  model  predicts  a  peak 
vedue  of  that  is  somewhat  lower  than  the  peak  value  of  =  4.5  suggested  by 
Patel  et  al.  (1985).  However,  this  is  a  known  characteristic  of  the  Launder- Sharma 
low-Reynolds-number  model,  tmd  is  of  little  consequence.  The  increase  in  turbulent 
kinetic  energy  at  the  interface  with  wind  speed  may  be  clearly  identified.  The 
profiles  are  similar  in  the  water,  except  that  the  increase  in  turbulent  kinetic  energy 
at  a  given  wind  speed  is  relatively  smaller,  since  the  flow  in  the  water  is  effectively 
smoother. 

Figure  7  shows  corresponding  mean  profiles  of  dissipation  rate.  The  characteris¬ 
tic  shape  of  the  plot  of  dissipation  rate  for  a  smooth  wall  is  evident,  and  the  profile 
agrees  extremely  well  with  the  measured  data  of  Laufer  (1954)  for  the  low  wind 
speed  (i.e.,  smooth  boundary)  czae.  For  higher  wind  speeds  the  value  of  dissipa¬ 
tion  rate  at  the  interface  increases  quickly,  then  begins  to  decrease  with  increasing 
wind  speed.  Also,  the  peak  in  the  dissipation  rate  at  about  y+  =  10  is  gradually 
reduced  as  wind  speed  increases.  The  model  correctly  predicts  a  non-zero  value  of 
c  at  the  interface,  in  contrast  to  the  low  turbulent  Reynolds  number  model  used  by 
Al-Zanaidi  and  Hui  (1984),  which  incorrectly  yields  e  =  0  at  a  boimdary. 
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5  Conclusions 


A  numerical  model  of  a  coupled,  stratified,  two-phase  flow  using  the  k-epsilon  tur¬ 
bulence  model  has  been  has  been  presented.  The  model  is  efficiently  solved  numer¬ 
ically  using  a  finite  difference  technique  with  Newton  iteration.  In  the  development 
of  the  model,  we  have  attempted  to  use  well  tested  methods  where  possible  (e.g.. 
the  k-epsilon  turbulence  model,  and  the  finite  difference  method  of  solution). 

The  numerical  model  has  been  applied  to  the  case  of  the  coupled  air-water  flow 
in  the  Stanford  Wind  Water  W'ave  Research  Facihty.  The  agreement  between  the 
mean  profiles  computed  by  numerical  model  and  experimental  data  is  encoiurag- 
ing.  However,  the  agreement  between  the  simulations  eind  the  experimental  data 
somewhat  better  in  the  air  than  in  the  water.  This  difference  reflects  our  neglect 
of  the  effects  of  wave- turbulence  interaction  in  the  water.  Regardless  of  this  defi¬ 
ciency,  the  model  predicts  the  behavior  of  the  mean  profiles  quantitatively,  as  well 
as  qualitatively.  In  addition,  the  method  for  treating  the  boimdary  conditions  at 
the  air-water  interface  shows  a  great  deal  of  promise. 

This  model  has  been  developed  as  part  of  a  larger  study  on  the  growth  of  wind¬ 
generated  waves.  The  mean  profiles  generated  by  the  numerical  model  will  be  used 
as  input  for  a  linear  stability  model  of  wave  growth  rate.  It  is  anticipated  that  the 
inclusion  of  both  the  shear  flow  in  the  water,  and  turbulent  effects  (via  the  k-epsilon 
turbulence  model)  in  this  model  of  wave  growth  will  result  in  improved  agreement 
between  predictions  and  meeisurements  of  wave  growth  rates,  especially  for  long 
waves. 


20 


References 


[1]  Akai,  M.,  A.  Inoue,  and  S.  Aoki,  “The  prediction  of  stratified  two-phase  flow 
with  a  two-equation  model  of  turbulence,”  Ini.  J.  Multiphase  Flow.,  Vol.  7,  pp. 
21-39,  1981. 

[2]  Al-Zanaidi,  M.A.,  and  W.H.  Hui,  “Turbulent  airflow  over  water  waves  -  a 
niunericaJ  study,”  J.  Fluid  Meek.,  Vol.  148,  pp.  225-246,  1984. 

[3]  Amorocho,  J.  and  J.J.  DeVries,  “A  new  evaluation  of  the  wind  stress  coefficient 
over  water  surfaces,”  J.  Geophys.  Research.,  Vol.  85(C1),  pp  433-442,  1980. 

[4]  Ascher,  U.M.,  R.M.M.  Mattheij,  and  R.D.  Russell,  “Numerical  solution  of 
boundary  value  problems  for  ordinary  differential  equations,”  Prentice  Hall, 
New  Jersey,  1988. 

[5]  Cheimg,  T.K.,  “A  study  of  the  turbulent  layer  in  the  water  at  an  air-water 
interface,”  Ph.D.  Dissertation  (Technictil  Report  No.  287),  Dept,  of  Civil  Eng., 
Stanford  Univ.,  1985. 

[6]  Cheung,  T.K.,  and  R.L.  Street,  “The  turbulent  layer  in  the  water  at  an  air- 
water  interface,”  J.  Fluid  Mech.,  Vol.  194,  pp.  133-151,  1988 

[7]  Coles,  D.,  “A  model  for  flow  in  the  viscous  sublayer,”  Proceedings  of  the  work¬ 
shop  on  coherent  structure  of  turbulent  boundary  layers.,  Lehigh  University, 
Bethlehem,  Pa.,  1978. 

[8]  Gent,  P.R.,  “A  numerical  model  of  the  air  flow  above  water  waves.  Part  2,”  J. 
Fluid  Mech.,  Vol.  82(2),  pp.  349-369,  1977. 

[9]  Hsu,  C.T.,  E.  Y.  Hsu,  and  R.L.  Street,  “On  the  structure  of  turbulent  flows 
over  a  progressive  water  wave:  theory  and  experiment  in  a  transformed,  wave¬ 
following  coordinate  system,”  J.  Fluid  Mech.,  Vol.  105,  pp.  87-117,  1981. 

[10]  Issa,  R.  I.,  “Prediction  of  turbulent,  stratified,  two-phase  flow  in  inclined  pipes 
and  channels,”  Ini.  J.  Multiphase  Flow.,  Vol.  14(2),  pp.  141-154,  1988. 

[11]  Kawai,  S.,  “Generation  of  initial  wavelets  by  instability  of  a  coupled  shear  flow 
and  their  evolution  to  wind  waves,”  J.  Fluid  Mech.,  Vol.  93(4),  pp.  661-703, 
1979. 

[12]  Launder,  B.E.,  and  B.I.  Sharma,  “Appheation  of  the  energy  dissipation  model 
of  turbulence  to  the  calculation  of  a  flow  near  a  spinning  disk,”  Letters  in  Heat 
and  Mass  Transfer.,  Vol.  1,  pp.  131-138,  1974. 

[13]  McIntosh,  D.A.,  R.L.  Street,  and  E.Y.  Hsu, “Turbulent  heat  and  momentum 
transfer  at  an  air- water  interfaK:e:  the  influence  of  surface  conditions,”  Technical 
Report  No.  197,  Dept,  of  Civil  Eng.,  Stanford  Univ.,  1975. 


21 


[14]  Miles,  J.VV.,  “On  the  generation  of  surfcice  waves  by  shear  flows."  J.  Fluid 
Mtch.,  Vol.  3,  pp.  185-204,  1957. 

[15]  Miles,  J.W.,  “On  the  generation  of  surface  waves  by  shear  flows,"  J.  Fluid 
Mech.,  Vol.  6,  pp.  469-478,  1959. 

[16]  Norris,  H.L.,  and  W.C.  Reynolds,  “Tmbulent  channel  flow  with  a  moving  wavy 
boundary,”  Technical  Report  TF-7,  Dept,  of  Mech.  Eng.,  Stanford  Univ.,  1975. 

[17]  Patankar,  S.V.,  and  D.B.  Spalding.,  “Heat  and  mass  transfer  in  boundary 
layers,”  Morgan-Grampian,  London,  1967. 

[18]  Patel,  V'.C.,  W.  Rodi,  and  G.  Scheuerer,  “Turbulence  models  for  near- wall 
and  low  Reynolds  number  flows:  A  review,”  AIAA  Journal.,  Vol.  23(9),  pp. 
1308-1319,  1985. 

[19]  Valenzuela,  G.R.,  “The  growth  of  gravity-capillary  waves  in  a  coupled  shear 
flow,”  J.  Fluid  Mech.,  Vol.  76(2),  pp.  229-250,  1976. 

[20]  van  Gastel,  K.,  P.A.E.M.  Janssen,  and  G.J.  Komen,  “On  phase  velocity  and 
growth  rate  of  wind-induced  gravity-ca^  U.  wave.s.”  J.  Fluid  Mech.,  Vol.  161, 
pp.  199-216,  1985. 

[21]  Zagustin,  K,  E.Y.  Hsu,  and  R.L  Street,  “Turbulent  flow  over  a  moving  bound¬ 
ary,”  Journal  of  the  Wateruji  ys  and  Harbors  Division,  ASCE,  Vol.  94,  No. 
WW4,  Proc.  Paper  6210,  pp.  397-414,  1968. 


22 


n/(  n  -  n) 


uAi 

•• 


Figure  3:  Mean  velocity  profile  in  water  for  Uoo  =  3.2  m/s 


Figure  4:  Mean  velocity  profiles  in  air  above  interface.  Comparison  of  numerical 
model  with  experiments  of  McIntosh  et  al.  (1975)  data,  o  -  2.3  m/s;  □  -  3.3  m/s;  o 
-  5.4  m/s;  x  -  8.1  m/s;  +  -  10.7  m/s 
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Figure  5:  Mean  velocity  profiles  in  air  above  interface.  Comparison  of  numerical 
model  with  experiments  of  McIntosh  et  al.  (1975)  data,  o  -  1.5  m/s;  □  -  2.6  m/s;  o 
-  3.2  m/s;  x  -  4.7  m/s;  +  -  6.7  m/s;  A  -  9.9  m/s;  •  -  13.1  m/s 
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